
# load required packages
pacman::p_load(ggplot2, tidyr, dplyr, readr, jtools, 
               ggeffects, fixest, interactions, texreg, 
               ipw, survey, zoo, ggpubr)

# Download data 
load("jcr_data.RData")

### Create conditioning sets 

# base controls
basecontrols <- c("lnyrsinceonst", "regionpronset",
                  "logoil", "logpop", "ldr_rebel",
                  "ldr_rchgcoup") 

# LDV + LDV X Year means + LDV X Unit means
LDV <- c("lallonset1_5_regm1", "lallonset1_5_yrm1",
         "mt1ldv", "lallonset1_5") 

# base controls regime-level means
baseunitm <- c("lnyrsinceonst_regm1",
                "regionpronset_regm1", "logoil_regm1", 
                "logpop_regm1", "ldr_rebel_regm1", 
                "ldr_rchgcoup_regm1") 

# base controls year-level means 
baseunityr <- c("lnyrsinceonst_yrm1", 
                 "regionpronset_yrm1", "logoil_yrm1",
                 "logpop_yrm1", "ldr_rebel_yrm1",
                 "ldr_rchgcoup_yrm1")

# unit-means X year-means of base controls + opposition in cabinet dummy 
unitXtime <- c("mt1", "mt2", "mt4",
                "mt5", "mt7", "mt8")

# unit and year means for opposition in cabinet dummy 
unitmyr1 <- c("lmultidum_regm1", "lmultidum_yrm1", "indvar1") 

# unit and year means for number of opposition groups
unitmyr2 <- c("logpartynr_1_regm1", "logpartynr_1_yrm1", "indvar2") 

# unit and year means for proportion of opposition elites
unitmyr3 <- c("logpercentopp_regm1", "logpercentopp_yrm1", "indvar3")

# unit and year means for weighted share of opposition elites
unitmyr4 <- c("logpercentopp_weighted_regm1", "logpercentopp_weighted_yrm1", "indvar4")

## conditioning sets with additional controls 

fullcontrols <- c("lnldryrsin", "lrgdpopct1", 
                  "lpolyarchy", "lpsinst",
                  "oppgain", "prior_civilwar") ## all controls 

baseunitmfull <- c("lnldryrsin_regm1", "lrgdpopct1_regm1", 
                    "lpolyarchy_regm1", "lpsinst_regm1",
                    "oppgain_regm1", "prior_civilwar_regm1") ## unit means (additional controls)

baseunityrfull <- c("lnyrsinceonst_yrm1",
                     "regionpronset_yrm1","lnldryrsin_yrm1", 
                     "lrgdpopct1_yrm1", "logoil_yrm1",
                     "logpop_yrm1", "ldr_rebel_yrm1", 
                     "ldr_rchgcoup_yrm1", "lpolyarchy_yrm1",
                     "lpsinst_yrm1", "oppgain_yrm1", 
                     "prior_civilwar_yrm1") ## year means (additional controls)

unitXtimefull <- c("mt3", "mt6", "mt9", "mt10", "mt11", 
                    "mt12") ## unit-means X year-means of additional controls 

### Opposition in Cabinet as IV

### Model 1 - Baseline 

jcr_data$lmultidum <- as.factor(jcr_data$lmultidum)

# Fit the model
mod1 <- glm(reformulate(c("lmultidum", basecontrols, baseunitm, baseunityr, unitmyr1, unitXtime), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data)

# Get clustered standard errors
mod1_clustered <- summ(mod1, cluster = "gwf_caseid", robust = "HC1", data = jcr_data)
print(mod1_clustered)

# Effect plot 
effectplot1 <- effect_plot(
  mod1, pred = "lmultidum", cluster = "gwf_caseid", robust = "HC1",
  int.width = 0.95, y.label = "Pr(Mass Uprising Onset)", x.label = "Opposition in Cabinet",
  interval = TRUE, cat.pred.point.size = 3, line.thickness = 2, point.shape = TRUE,
  cat.geom = "point", cat.interval.geom = "linerange") + 
  ylim(0, 0.1) +
  theme_classic() +
  theme(axis.text.x = element_text(color = "black", size = 25), 
    axis.text.y = element_text(color = "black", size = 25),
    axis.title.x = element_text(size = 25), 
    axis.title.y = element_text(size = 25))

effectplot1

### Model 2 - LDV (Dummy for uprising onset in the past 5 year)
mod1.ldv <- glm(reformulate(c("lmultidum", basecontrols, baseunitm, baseunityr, unitmyr1, unitXtime, LDV),
                            response = "allonset"), family = binomial(link = "probit"), data = jcr_data)

mod1.ldv <- summ(mod1.ldv, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod1.ldv

### Model 3 - Additional controls 
mod1.full1 <- glm(reformulate(c("lmultidum", basecontrols, baseunitm, baseunityr, unitmyr1, unitXtime, LDV, 
                                fullcontrols, baseunitmfull, baseunityrfull, unitXtimefull), 
                              response = "allonset"), family = binomial(link = "probit"), data = jcr_data)

mod1.full1 <- summ(mod1.full1, cluster = "gwf_caseid", robust = "HC1", data = jcr_data)
mod1.full1

### Number of Opposition Groups in Cabinet as IV 

### Model 5 - Baseline 
mod2 <- glm(reformulate(c("log(lpartynr + 1)", basecontrols, baseunitm, baseunityr, unitmyr2, unitXtime), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod2_clustered <- summ(mod2, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod2_clustered

effectplot2 <- effect_plot(mod2, pred = lpartynr, cluster = "gwf_caseid", robust = "HC1", 
                           data = jcr_data, interval = TRUE, rug = TRUE, 
                           rug.sides = "b", jitter = 0.3, line.thickness = 1, 
                           y.label = "Pr(Mass Uprising Onset)") + 
  scale_x_continuous(name = "Number of Opposition Groups in Cabinet", 
                     breaks = 0:7, limits = c(0, 7)) +
  ylim(0, 0.1) + 
  theme_classic() +
  theme(axis.text.x = element_text(color = "black", size = 25), 
    axis.text.y = element_text(color = "black", size = 25),
    axis.title.x = element_text(size = 25),
    axis.title.y = element_text(size = 25))

effectplot2

### Model 6 - LDV 

mod2.ldv <- glm(reformulate(c("log(lpartynr + 1)", basecontrols, baseunitm, baseunityr, unitmyr2, unitXtime, LDV), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 
mod2.ldv <- summ(mod2.ldv, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod2.ldv

### Model 7 - Additional controls 

mod2.full2 <- glm(reformulate(c("log(lpartynr + 1)", basecontrols, baseunitm, baseunityr, unitmyr2, unitXtime, LDV, 
                              fullcontrols, baseunitmfull, baseunityrfull, unitXtimefull), 
                            response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 
mod2.full2 <- summ(mod2.full2, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod2.full2

### Proportion of Opposition Elites in Cabinet as IV 

### Model 9 - Baseline 
mod3 <- glm(reformulate(c("log((lprop_opp_member*100) + 1)", basecontrols, baseunitm, baseunityr, unitmyr3, unitXtime), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod3_clustered <- summ(mod3, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod3_clustered

effectplot3 <- effect_plot(mod3, pred = lprop_opp_member, cluster = "gwf_caseid", robust = "HC1", 
                           data = jcr_data, interval = TRUE, rug = TRUE, rug.sides = "b", 
                           line.thickness = 1, y.label = "Pr(Mass Uprising Onset)", 
                           x.label = "Proportion of Opposition Elites in Cabinet") + 
  ylim(0, 0.1) + 
  theme_classic() +
  theme(axis.text.x = element_text(color = "black", size = 25),
        axis.text.y = element_text(color = "black", size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))

effectplot3

### Model 10 - LDV 

mod3.ldv <- glm(reformulate(c("log((lprop_opp_member*100) + 1)", basecontrols, baseunitm, baseunityr, unitmyr3, unitXtime, LDV), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod3.ldv <- summ(mod3.ldv, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod3.ldv

### Model 11 - Additional controls 

mod3.full3 <- glm(reformulate(c("log((lprop_opp_member*100) + 1)", basecontrols, baseunitm, baseunityr, unitmyr3, unitXtime, LDV, 
                                fullcontrols, baseunitmfull, baseunityrfull, unitXtimefull), 
                                response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod3.full3 <- summ(mod3.full3, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod3.full3

### Weighted Share of Opposition Elites in Cabinet as IV 

### Model 13 - Baseline 
mod4 <- glm(reformulate(c("log((lweighted_opposition_share*100) + 1)", basecontrols, baseunitm, baseunityr, unitmyr4, unitXtime), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod4_clustered <- summ(mod4, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod4_clustered

effectplot4 <- effect_plot(mod4, pred = lweighted_opposition_share, cluster = "gwf_caseid", robust = "HC1", 
                           data = jcr_data, interval = TRUE, rug = TRUE, rug.sides = "b", 
                           line.thickness = 1, y.label = "Pr(Mass Uprising Onset)", 
                           x.label = "Proportion of Opposition Elites in Cabinet - Weighted") + 
  ylim(0, 0.1) + 
  theme_classic() +
  theme(axis.text.x = element_text(color = "black", size = 25),
    axis.text.y = element_text(color = "black", size = 25),
    axis.title.x = element_text(size = 25),
    axis.title.y = element_text(size = 25))

effectplot4

### Model 14 - LDV

mod4.ldv <- glm(reformulate(c("log((lweighted_opposition_share*100) + 1)", 
                              basecontrols, baseunitm, baseunityr, unitmyr4, unitXtime, LDV), 
                        response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod4.ldv <- summ(mod4.ldv, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod4.ldv

### Model 15 - All controls 

mod4.full4 <- glm(reformulate(c("log((lweighted_opposition_share*100) + 1)", 
                                basecontrols, baseunitm, baseunityr, unitmyr4, unitXtime, LDV, 
                                fullcontrols, baseunitmfull, baseunityrfull, unitXtimefull), 
                                response = "allonset"), family = binomial(link = "probit"), data = jcr_data) 

mod4.full4 <- summ(mod4.full4, cluster = "gwf_caseid", robust = "HC1", data = jcr_data) 
mod4.full4

### Marginal Structural Models 

jcr_data$lmultidum <- as.numeric(jcr_data$lmultidum)

# Data preparation

# Create dataset for IPW: ipwtm needs data frame without missing values
foripw_1 <- jcr_data %>% 
  dplyr::select(
  gwf_caseid, year, lnyrsinceonst, regionpronset, lnldryrsin, lrgdpopct1, logoil, logpartynr_1, logpartynr_1_regm1, logpartynr_1_yrm1,
                                             logpercentopp_regm1, logpercentopp_yrm1, logweightedshare_regm1, logweightedshare_yrm1, alesinaeth, alesinaeth_dum_med, alesinaeth_dum_q3, fearoneth_dum_q3, fearoneth_dum_med, fearoneth, 
                                             logpop, lpolyarchy, regionid, lpsinst, lmultidum, laglogpercentopp, laglogweightedshare, laglogpartynr1, laglmultidum1, lmultidum1_5, 
                                              logweightedshare, allonset1_5, lmultidum_regm1, lmultidum_yrm1, v2psoppaut, v2psparban, lv2xpe_exlsocgr, collegialparty,
                                              oppgain, logpercentopp, allonset, ldr_rebel, ldr_rchgcoup, pcoldwar,
                                              lpsinst, prior_civilwar) %>% 
    drop_na()

table(foripw_1$lmultidum)

# Weight calculation using IPW
w.dum.ipw1 <- ipwtm(
  exposure = lmultidum,
  family = "binomial",
  link = "logit",
  numerator = ~ laglmultidum1 + lmultidum1_5 + alesinaeth_dum_q3 + ldr_rebel + ldr_rchgcoup + prior_civilwar + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglmultidum1 + lmultidum1_5 + allonset1_5  + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst +   
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar),
  id = gwf_caseid,
  timevar = year,
  type = "all",
  data = as.data.frame(foripw_1))

ipwplot(weights = w.dum.ipw1$ipw.weights, timevar = foripw_1$year, binwidth=1, 
                   main="", ylab="(Log)Mean of Weights", xlab="Consecutive Observations (Years)")

summary(w.dum.ipw1$ipw.weights)

foripw_1$sw1 <- w.dum.ipw1$ipw.weights

### Model 4 - MSM
msm_1 <- svyglm(allonset ~ lmultidum + lmultidum_regm1 + lmultidum_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw1, 
                                   data = foripw_1), family = "binomial")
summary(msm_1)

### Model 8 - MSM 

w.pnr.ipw2 <- ipwtm(
  exposure = logpartynr_1,
  family = "gaussian",
  numerator = ~ laglogpartynr1 + lmultidum1_5 + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglogpartynr1 + lmultidum1_5 + allonset1_5  + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

foripw_1$sw2 <- w.pnr.ipw2$ipw.weights

summary(w.pnr.ipw2$ipw.weights)

ipwplot(weights = w.pnr.ipw2$ipw.weights, timevar = foripw_1$year, binwidth=1, 
                   main="", ylab="(Log)Mean of Weights", xlab="Consecutive Observations (Years)")

msm_2 <- svyglm(formula = allonset ~ logpartynr_1 + logpartynr_1_regm1 + logpartynr_1_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw2, 
                                   data = foripw_1), family = "binomial")
summary(msm_2)

### Model 12 - MSM 

w.pct.ipw <- ipwtm(
  exposure = logpercentopp,
  family = "gaussian",
  numerator = ~ laglogpercentopp + lmultidum1_5 + alesinaeth_dum_q3 + ldr_rebel + ldr_rchgcoup + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  denominator = ~ laglogpercentopp + lmultidum1_5 + allonset1_5 + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

summary(w.pct.ipw$ipw.weights)

ipwplot(weights = w.pct.ipw$ipw.weights, timevar = foripw_1$year, binwidth=1, 
                   main="", ylab="(Log)Mean of Weights", xlab="Consecutive Observations (Years)")

foripw_1$sw3 <- w.pct.ipw$ipw.weights

msm_3 <- svyglm(allonset ~ logpercentopp + logpercentopp_regm1 + logpercentopp_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw3, 
                                   data = foripw_1), family = "binomial")
summary(msm_3)

### Model 16 - MSM 

w.pct.ipw_weighted <- ipwtm(
  exposure = logweightedshare,
  family = "gaussian",
  numerator = ~ laglogweightedshare + lmultidum1_5 + alesinaeth_dum_q3 + ldr_rebel + ldr_rchgcoup + prior_civilwar + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglogweightedshare + lmultidum1_5 + allonset1_5 + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

summary(w.pct.ipw_weighted$ipw.weights)

foripw_1$sw4 <- w.pct.ipw_weighted$ipw.weights

ipwfig4 <- ipwplot(weights =  w.pct.ipw_weighted$ipw.weights, timevar = foripw_1$year, binwidth=1, 
                   main="", ylab="(Log)Mean of Weights", xlab="Consecutive Observations (Years)")

msm_4 <- svyglm(allonset ~ logweightedshare + logweightedshare_regm1 + logweightedshare_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw4, 
                                   data = foripw_1), family = "binomial")
summary(msm_4)

## Coefficient Plot (Figure 2)
model_list <- list(mod1_clustered, mod1.ldv, mod1.full1, msm_1, 
                   mod2_clustered, mod2.ldv, mod2.full2, msm_2, 
                   mod3_clustered, mod3.ldv, mod3.full3, msm_3, 
                   mod4_clustered, mod4.ldv, mod4.full4, msm_4)

coefplot1 <- plot_summs(model_list, 
                        point.shape = c(0, 1, 2, 5, 6, 7, 8, 9, 3, 4, 10 ,13, 15, 16 , 17 , 18), 
                        ci_level = 0.95, inner_ci_level = 0.90, 
                        colors = viridis(length(model_list)), 
                        point.size = 10, line.size = c(3, 4), 
                        model.names = c("Model 1 - Baseline", "Model 2 - LDV ","Model 3 - All controls",  "Model 4 - MSM",
                                        "Model 5 - Baseline ", "Model 6 - LDV ", "Model 7 - All controls", "Model 8 - MSM", 
                                        "Model 9 - Baseline", "Model 10 - LDV", "Model 11 - All controls", "Model 12 - MSM", 
                                        "Model 13 - Baseline", "Model 14 - LDV", "Model 15 - All controls", "Model 16 - MSM"), 
                        coefs = c("Opposition (Dummy)"  = "lmultidum", 
                                "Ln(No. of Opposition Groups)" = "log(lpartynr + 1)",
                                "Ln(No. of Opposition Groups)" = "logpartynr_1", 
                                "Ln(Prop. Opposition Elites)" = "logpercentopp", 
                                "Ln(Prop. Opposition Elites)" = "log((lprop_opp_member * 100) + 1)", 
                                "Ln(Prop. Opposition Elites - Weighted)" = "log((lweighted_opposition_share * 100) + 1)",
                                "Ln(Prop. Opposition Elites - Weighted)" = "logweightedshare"), exp = F) + 
                        labs(title = "All Uprisings", x = "Coefficient Estimates") + 
                        theme_classic()  


coefplot1 <- coefplot1  + theme(plot.title = element_blank(),
                                legend.text = element_text(size=30), 
                                legend.position="right", 
                                legend.title = element_blank(), 
                                axis.title.y = element_blank(),
                                axis.title.x = element_text(size = 32),
                                axis.text.x = element_text(color="black", size=32), 
                                axis.text.y = element_text(color="black", size=32), 
                                legend.margin = margin(0, 0, 0, 0),
                                legend.spacing.x = unit(17, "mm"), 
                                legend.key.height = unit(1, 'cm'), 
                                legend.key.width = unit(3, 'cm'),
                                legend.spacing.y = unit(30, "mm")) + 
                                guides(color=guide_legend(ncol=1, nrow=16)) + 
                                scale_x_continuous(breaks = c(-2, -1.5, -1, -0.5, 0),
                                                   limits = c(-2, 0))

ggsave("coefplot1.pdf", coefplot1, width=25, height=10, units="in")


# Figure 3 
effectplotsmerged <- ggarrange(effectplot1, effectplot2, effectplot3, effectplot4, ncol = 1, nrow = 4, 
                               labels = c("A", "B", "C", "D"), font.label=list(color="black",size=22))
effectplotsmerged
ggsave("effectplotsmerged.pdf", effectplotsmerged, width=10, height=25, units="in")


## Main regression results in table
texreg(list(mod1_clustered, mod1.ldv, mod1.full1, msm_1),stars = c(0.01,0.05, 0.1))
texreg(list(mod2_clustered, mod2.ldv, mod2.full2, msm_2),stars = c(0.01,0.05, 0.1))
texreg(list(mod3_clustered, mod3.ldv, mod3.full3, msm_3),stars = c(0.01,0.05, 0.1))
texreg(list(mod4_clustered, mod4.ldv, mod4.full4, msm_4),stars = c(0.01,0.05, 0.1))

# Line plots (Figure 1)

jcr_data$lmultidum <- as.numeric(prcabdata$jcr_data)


lineplot <- jcr_data %>% group_by(year) %>% mutate(globalpropopp = mean(lprop_opp_member), 
                                                   globalpropopp_w = mean(lweighted_opposition_share),
                                                   globalpartynr = mean(lpartynr),
                                                   unique_cc = n_distinct(ccode), 
                                                   countmulti = sum(lmultidum), 
                                                   globalmultidum = countmulti/unique_cc) %>% 
                                                   dplyr::select(globalpropopp, globalpropopp_w, countmulti, globalpartynr, globalmultidum,
                                                                   unique_cc, year) %>% distinct() 

base_theme <- theme_classic() +
  theme(
    plot.title = element_text(size = 11),
    axis.text.x = element_text(color = "black", size = 12),
    axis.text.y = element_text(color = "black", size = 12),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

# Helper function to create line plots
create_lineplot <- function(data, y_var, title, y_limits, y_breaks, y_label = NULL) {
  ggplot(data, aes(x = year, y = !!sym(y_var))) +
    geom_line(alpha = 0.3) +
    geom_smooth(color = "black", span = 1) +
    scale_y_continuous(name = y_label, limits = y_limits, breaks = y_breaks) +
    scale_x_continuous(name = "", breaks = c(1970, 1980, 1990, 2000, 2010, 2020)) +
    labs(title = title) +
    base_theme
}

lineplot1 <- create_lineplot(lineplot, "globalmultidum", "% of regimes with the opposition in cabinet", 
                             c(0, 0.8), c(0, 0.2, 0.4, 0.6, 0.8))

lineplot2 <- create_lineplot(lineplot, "globalpartynr", "Mean number of opposition groups in cabinets", 
                             c(0, 3), c(0, 1, 2, 3), "Mean opposition groups in cabinets")

lineplot3 <- create_lineplot(lineplot, "globalpropopp", "% of opposition elites in cabinets", 
                             c(0, 0.3), c(0, 0.1, 0.2, 0.3), "Mean proportion of opposition elites in cabinets")

lineplot4 <- create_lineplot(lineplot, "globalpropopp_w", "% of opposition elites in cabinets - weighted", 
                             c(0, 0.3), c(0, 0.1, 0.2, 0.3))

lineplots_merged <- ggarrange(
  lineplot1, lineplot2, lineplot3, lineplot4,
  ncol = 2, 
  labels = c("A", "B", "C", "D"),
  font.label = list(size = 10),
  nrow = 2)

lineplots_merged

ggsave("lineplotssmerged2.pdf", lineplots_merged, width=8, height=5, units="in")

## Robustness checks 

# Sensitivity Analyses (Section 6.1 Appendix)

# Opposition in cabinet dummy as IV
mod_delta1 <- o_delta(y = "allonset", x = "lmultidum", 
                      con = "lnyrsinceonst + regionpronset + 
logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", 
                      id = "gwf_caseid", time = "year", beta = 0,
                      R2max = 0.153 * 1.3, type = "lm", data = jcr_data)
mod_delta1

mod_delta2 <- o_delta(y = "allonset", x = "lmultidum", 
                      con = "lnyrsinceonst + regionpronset + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
    oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", 
                      id = "gwf_caseid", time = "year", beta = 0,
                      R2max = 0.153*1.7, type = "lm", data = jcr_data)
mod_delta2

mod_delta3 <- o_delta(y = "allonset", x = "lmultidum", 
                      con = "lnyrsinceonst + regionpronset + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
    oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", 
                      id = "gwf_caseid", time = "year", beta = 0,
                      R2max = 0.153*2, type = "lm", data = jcr_data)

mod_delta3

mod_delta4 <- o_delta(y = "allonset", x = "lmultidum", 
                      con = "lnyrsinceonst + regionpronset + 
logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
 oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", 
                      id = "gwf_caseid", time = "year", beta = 0,
                      R2max = 0.153*2.6, type = "lm", data = jcr_data)
mod_delta4

sensplot1 <- o_delta_rsq_viz(y = "allonset",           
                              x = "lmultidum", con = "lnyrsinceonst + regionpronset + 
                  logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                  oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5",
                              beta = 0,                    
                              type = "lm",               
                              data = jcr_data)    

sensplot1 <- sensplot1 + 
  scale_y_continuous(breaks = 0:10, limits = c(0, 10)) + 
  scale_x_continuous(breaks = seq(0.2, 0.5, by = 0.1), limits = c(0.15, 0.5)) + 
  labs(title = "Sensitivity Test for Opposition Group (Dummy)") + 
  theme_classic() + 
  theme(
    plot.title = element_text(size = 17),
    axis.text = element_text(color = "black", size = 15), 
    axis.title = element_text(size = 17))
sensplot1

# Number of opposition groups as IV

mod2_delta1 <- o_delta(y = "allonset", x = "logpartynr_1", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.3, type = "lm", data = jcr_data)
mod2_delta1

mod2_delta2 <- o_delta(y = "allonset", x = "logpartynr_1", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.7, type = "lm", data = jcr_data)
mod2_delta2
mod2_delta3 <- o_delta(y = "allonset", x = "logpartynr_1", con = "lnyrsinceonst + regionpronset + 
                         logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                         oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2, type = "lm", data = jcr_data)
mod2_delta3
mod2_delta4 <- o_delta(y = "allonset", x = "logpartynr_1", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2.6, type = "lm", data = jcr_data)
mod2_delta4

sensplot2 <- o_delta_rsq_viz(y = "allonset",           
                               x = "logpartynr_1", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5",
                               beta = 0,                    
                               type = "lm",               
                               data = jcr_data)  

sensplot2 <- sensplot2 + 
  scale_y_continuous(breaks = 0:10, limits = c(0, 10)) + 
  scale_x_continuous(breaks = seq(0.2, 0.5, by = 0.1), limits = c(0.15, 0.5)) + 
  labs(title = "Sensitivity Test for Ln(N. Opposition Groups)") + 
  theme_classic() + 
  theme(
    plot.title = element_text(size = 17),
    axis.text = element_text(color = "black", size = 15), 
    axis.title = element_text(size = 17))

sensplot2

# Proportion of opposition elites as IV

mod3_delta1 <- o_delta(y = "allonset", x = "logpercentopp", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.3, type = "lm", data = jcr_data)
mod3_delta1

mod3_delta2 <- o_delta(y = "allonset", x = "logpercentopp", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.7, type = "lm", data = jcr_data)
mod3_delta2

mod3_delta3 <- o_delta(y = "allonset",  x = "logpercentopp",  con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2, type = "lm", data = jcr_data)
mod3_delta3

mod3_delta4 <- o_delta(y = "allonset",  x = "logpercentopp",  con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2.6, type = "lm", data = jcr_data)
mod3_delta4


sensplot3 <- o_delta_rsq_viz(y = "allonset",           
                               x = "logpercentopp", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5",
                               beta = 0,                    
                               type = "lm",               
                               data = jcr_data)  

sensplot3 <- sensplot3 + 
  scale_y_continuous(breaks = 0:10, limits = c(0, 10)) + 
  scale_x_continuous(breaks = seq(0.2, 0.5, by = 0.1), limits = c(0.15, 0.5)) + 
  labs(title = "Sensitivity Test for Ln(Prop. Opposition Elites)") + 
  theme_classic() + 
  theme(
    plot.title = element_text(size = 17),
    axis.text = element_text(color = "black", size = 15), 
    axis.title = element_text(size = 17))

sensplot3

# Weighted proportion of opposition elites as IV

mod4_delta1 <- o_delta(y = "allonset", x = "logweightedshare", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.3, type = "lm", data = jcr_data)
mod4_delta1
mod4_delta2 <- o_delta(y = "allonset", x = "logweightedshare", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*1.7, type = "lm", data = jcr_data)
mod4_delta2
mod4_delta3 <- o_delta(y = "allonset",  x = "logweightedshare", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2, type = "lm", data = jcr_data)
mod4_delta3
mod4_delta4 <- o_delta(y = "allonset",  x = "logweightedshare", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5", id = "gwf_caseid", time = "year", beta = 0,
                       R2max = 0.153*2.6, type = "lm", data = jcr_data)
mod4_delta4

sensplot4 <- o_delta_rsq_viz(y = "allonset",           
                               x = "logweightedshare", con = "lnyrsinceonst + regionpronset + 
                        logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpsinst + 
                        oppgain +  factor(gwf_caseid) + factor(year) + lallonset1_5",
                               beta = 0,                    
                               type = "lm",               
                               data = jcr_data)  

sensplot4 <- sensplot4 + 
  scale_y_continuous(breaks = 0:10, limits = c(0, 10)) + 
  scale_x_continuous(breaks = seq(0.2, 0.5, by = 0.1), limits = c(0.15, 0.5)) + 
  labs(title = "Sensitivity Test for Ln(Prop. Opposition Elites - Weighted)") + 
  theme_classic() + 
  theme(
    plot.title = element_text(size = 17),
    axis.text = element_text(color = "black", size = 15), 
    axis.title = element_text(size = 17))

sensplot4

combinesens <- ggarrange(sensplot1, sensplot2, sensplot3, sensplot4, ncol = 2, nrow = 2, 
                         labels = c("A", "B", "C", "D"))
combinesens

ggsave("combinesens.pdf",combinesens, width=20, height=7, units="in")

## Testing for pre-trends (Table 8 Appendix)

variables <- c("lmultidum", "logpartynr_1", "logpercentopp", "logweightedshare")
leads <- c(1, 3, 5) 

models <- list()

for (var in variables) {
  for (lead in leads) {
    # Define the formula dynamically for leads
    formula <- as.formula(paste("allonset ~ f(", var, ", ", lead, ":", lead, ") + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid + year", sep = ""))
    
    # Fit the model and store it in the list
    mod <- feglm(formula, panel.id = ~gwf_caseid + year, data = jcr_data, family = binomial(link = "logit"))
    models[[paste(var, lead, sep = "_")]] <- mod
  }
}


# Use texreg to display the results
texreg(models, stars = c(0.01, 0.05, 0.1))

# Balance check (Appendix Section 6.3)

foripw_1$lmultidum <- as.numeric(foripw_1$lmultidum)


check_bal1 <- svyglm(formula = lmultidum ~ laglmultidum1 + lmultidum1_5 + allonset1_5  + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst +  
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, weights = ~sw1, data = foripw_1)) 
summary(check_bal1)

check_bal2 <- svyglm(formula = lmultidum ~ laglmultidum1 + lmultidum1_5 + allonset1_5 + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, data = foripw_1)) 
summary(check_bal2)


check_bal3 <- svyglm(formula = logpartynr_1 ~ laglogpartynr1 + lmultidum1_5 + allonset1_5  + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst +  
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, weights = ~sw2, data = foripw_1)) 

summary(check_bal3)

check_bal4 <- svyglm(formula = logpartynr_1 ~ laglogpartynr1 + lmultidum1_5 + allonset1_5 + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, data = foripw_1)) 
summary(check_bal4)

check_bal5 <- svyglm(formula = logpercentopp ~ laglogpercentopp  + lmultidum1_5 + allonset1_5  +
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst +  
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, weights = ~sw3, data = foripw_1)) 
summary(check_bal5)

check_bal6 <- svyglm(formula = logpercentopp ~ laglogpercentopp + lmultidum1_5 + allonset1_5 + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, data = foripw_1)) 

summary(check_bal6)

check_bal7 <- svyglm(formula = logweightedshare ~ laglogweightedshare + lmultidum1_5 + allonset1_5 + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, weights = ~sw4, data = foripw_1)) 

summary(check_bal7)

check_bal8 <- svyglm(formula = logweightedshare ~ laglogweightedshare + lmultidum1_5 + allonset1_5 + 
                       regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
                       logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + 
                       factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
                     design = svydesign(~gwf_caseid, data = foripw_1)) 

summary(check_bal8)

# Table 9 Appendix 
texreg(list(check_bal2, check_bal1, check_bal4, check_bal3, 
            check_bal6, check_bal5, check_bal8, check_bal7), 
      stars = c(0.01,0.05, 0.1))

# Alternative Specification for Weights  

w.dum.ipw1_alt <- ipwtm(
  exposure = lmultidum,
  family = "binomial",
  link = "logit",
  numerator = ~ laglmultidum1 + lmultidum1_5  + alesinaeth_dum_q3 + ldr_rebel + ldr_rchgcoup + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglmultidum1 + lmultidum1_5 + allonset1_5 + regionpronset + ldr_rebel + ldr_rchgcoup  + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + lpolyarchy + lpolyarchy + lpsinst + v2psoppaut + v2psparban +  lv2xpe_exlsocgr + collegialparty + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  data = as.data.frame(foripw_1))

foripw_1$sw1_alt <- w.dum.ipw1_alt$ipw.weights

msm_1_alt <- svyglm(allonset ~ lmultidum + lmultidum_regm1 + lmultidum_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw1_alt, 
                                   data = foripw_1), family = "binomial")
summary(msm_1_alt)

w.pnr.ipw2_alt <- ipwtm(
  exposure = logpartynr_1,
  family = "gaussian",
  numerator = ~ laglogpartynr1 + lmultidum1_5 + alesinaeth_dum_q3 + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglogpartynr1 + lmultidum1_5 + allonset1_5  + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + v2psoppaut + v2psparban +  lv2xpe_exlsocgr + collegialparty + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

foripw_1$sw2_alt <- w.pnr.ipw2_alt$ipw.weights

msm_2_alt <- svyglm(formula = allonset ~ logpartynr_1 + logpartynr_1_regm1 + logpartynr_1_yrm1, 
                design = svydesign(~gwf_caseid, weights = ~sw2_alt, 
                                   data = foripw_1), family = "binomial")
summary(msm_2_alt)

w.pct.ipw_alt <- ipwtm(
  exposure = logpercentopp,
  family = "gaussian",
  numerator = ~ laglogpercentopp + lmultidum1_5 + alesinaeth_dum_q3 +prior_civilwar+ factor(regionid) + factor(pcoldwar), 
  denominator = ~ laglogpercentopp + lmultidum1_5 + allonset1_5 + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + v2psoppaut + v2psparban +  lv2xpe_exlsocgr + collegialparty + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

foripw_1$sw3_alt <- w.pct.ipw_alt$ipw.weights

msm_3_alt <- svyglm(formula = allonset ~ logpercentopp + logpercentopp_regm1 + logpercentopp_yrm1, 
                    design = svydesign(~gwf_caseid, weights = ~sw3_alt, 
                                       data = foripw_1), family = "binomial")
summary(msm_3_alt)

w.pct.ipw_weighted_alt <- ipwtm(
  exposure = logweightedshare,
  family = "gaussian",
  numerator = ~ laglogweightedshare + lmultidum1_5 + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar),
  denominator = ~ laglogweightedshare + lmultidum1_5 + allonset1_5 + regionpronset + lnyrsinceonst + lnldryrsin + lrgdpopct1 + 
    logoil + logpop + ldr_rebel + ldr_rchgcoup + lpolyarchy + lpolyarchy + lpsinst + v2psoppaut + v2psparban +  lv2xpe_exlsocgr + collegialparty + 
    factor(oppgain) + alesinaeth_dum_q3 + prior_civilwar + factor(regionid) + factor(pcoldwar), 
  id = gwf_caseid,
  timevar = year,
  type = "all",
  corstr = "ar1",
  data = as.data.frame(foripw_1))

foripw_1$sw4_alt <- w.pct.ipw_weighted_alt$ipw.weights

msm_4_alt <- svyglm(formula = allonset ~ logweightedshare + logweightedshare_regm1 + logweightedshare_yrm1, 
                    design = svydesign(~gwf_caseid, weights = ~sw4_alt, 
                                       data = foripw_1), family = "binomial")
summary(msm_4_alt)

# Table 10 Appendix 
texreg(list(msm_1_alt, msm_2_alt, msm_3_alt, msm_4_alt), stars = c(0.01,0.05, 0.1))

# Alternative estimators Appendix Section 6.4 

altest1 <- feols(allonset ~ lmultidum + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], data = jcr_data)
summary(altest1)

altest2 <- feglm(allonset ~ lmultidum + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                 data = jcr_data, 
                 family = binomial(link = "logit"))
summary(altest2)

altest3 <- feols(allonset ~ logpartynr_1 + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], data = jcr_data)
summary(altest3)

altest4 <- feglm(allonset ~ logpartynr_1 + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                 data = jcr_data, 
                 family = binomial(link = "logit"))
summary(altest4)

altest5 <- feols(allonset ~ logpercentopp + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], data = jcr_data)
summary(altest5)

altest6 <- feglm(allonset ~ logpercentopp + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                 data = jcr_data, 
                 family = binomial(link = "logit"))
summary(altest6)

altest7 <- feols(allonset ~ logweightedshare + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], data = jcr_data)
summary(altest7)

altest8 <- feglm(allonset ~ logweightedshare + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                 data = jcr_data, 
                 family = binomial(link = "logit"))
summary(altest8)

# Table 11 Appendix
texreg(list(altest2, altest1, altest4, altest3, altest6, altest5, altest8, altest7),stars = c(0.01,0.05, 0.1))

# Excluding postconflict settings section 6.5 Appendix 
mod1.post1 <- feglm(allonset ~ lmultidum + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                    data = prcabdata_rev[prcabdata_rev$postconflict != 1,], 
                    family = binomial(link = "logit"))
summary(mod1.post1)
mod2.post2 <- feglm(allonset ~ logpartynr_1 + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                    data = prcabdata_rev[prcabdata_rev$regionid != 1,], 
                    family = binomial(link = "logit"))
summary(mod2.post2)
mod3.post3 <- feglm(allonset ~ logpercentopp + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                    data = prcabdata_rev[prcabdata_rev$regionid != 1,], 
                    family = binomial(link = "logit"))
summary(mod3.post3)
mod4.post4 <- feglm(allonset ~ logweightedshare + lnyrsinceonst + regionpronset  + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year], 
                    data = prcabdata_rev[prcabdata_rev$regionid != 1,], 
                    family = binomial(link = "logit"))
summary(mod4.post4)

# Table 12 Appendix
texreg(list(mod1.post1, mod2.post2, mod3.post3, mod4.post4),stars = c(0.01,0.05, 0.1))

# Excluding each region section 6.6 Appendix 

regions_to_exclude <- 0:6

# Table 13 Appendix 

mod1_results <- lapply(regions_to_exclude, function(region) {
  subset_data <- jcr_data[jcr_data$regionid != region, ]
  feglm(fml = allonset ~ lmultidum + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year],
        data = subset_data,
        family = binomial(link = "logit"))})

# Table 13 Appendix 
texreg(mod1_results,
       stars = c(0.01, 0.05, 0.1),
       caption = "Regression Models Excluding Regions (logpercentopp)",
       custom.model.names = paste("Excluding Region", regions_to_exclude))

mod2_results <- lapply(regions_to_exclude, function(region) {
  subset_data <- jcr_data[jcr_data$regionid != region, ]
  feglm(fml = allonset ~ logpartynr_1 + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year],
        data = subset_data,
        family = binomial(link = "logit"))})

#Table 14 Appendix 
texreg(mod2_results,
       stars = c(0.01, 0.05, 0.1),
       caption = "Regression Models Excluding Regions (logpercentopp)",
       custom.model.names = paste("Excluding Region", regions_to_exclude))


mod3_results <- lapply(regions_to_exclude, function(region) {
  subset_data <- jcr_data[jcr_data$regionid != region, ]
    feglm(fml = allonset ~ logpercentopp + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year],
    data = subset_data,
    family = binomial(link = "logit"))})

#Table 15 Appendix 
texreg(mod3_results,
  stars = c(0.01, 0.05, 0.1),
  caption = "Regression Models Excluding Regions (logpercentopp)",
  custom.model.names = paste("Excluding Region", regions_to_exclude))

mod4_results <- lapply(regions_to_exclude, function(region) {
  subset_data <- jcr_data[jcr_data$regionid != region, ]
  feglm(fml = allonset ~ logweightedshare + lnyrsinceonst + regionpronset + logoil + logpop + ldr_rebel + ldr_rchgcoup | gwf_caseid[year],
    data = subset_data,
    family = binomial(link = "logit"))})


# Table 16 Appendix 

texreg(mod4_results,
  stars = c(0.01, 0.05, 0.1),
  caption = "Regression Models Excluding Regions",
  custom.model.names = paste("Excluding Region", regions_to_exclude))



